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Abstract.  The  evaporation  of  a  liquid  slab  into  vacuum  is  studied  by  numerical  solutions  of  the  Enskog- Vlasov  equation  for 
a  fluid  of  spherical  molecules  interacting  by  Sutherland  potential.  The  main  aim  of  this  work  is  to  obtain  the  structure  of  the 
vapor-liquid  interface  in  non-equilibrium  conditions  as  well  as  the  distribution  function  of  evaporating  molecules.  The  results 
show  that  the  distribution  function  of  molecules  crossing  a  properly  defined  vapor-liquid  boundary  is  almost  Maxwellian  and 
that  the  vapor  phase  is  reasonably  well  described  by  the  Boltzmann  equation  with  diffusive  boundary  condition. 


INTRODUCTION 


In  many  situations  of  great  theoretical  and  practical  interest  it  is  necessary  to  study  the  motion  of  a  vapor  in  contact 
with  its  condensed  phase.  Most  studies  of  evaporation/condensation  phenomena  adopt  kinetic  theory  methods  and 
focus  on  the  behavior  of  the  vapor  phase[l].  Accordingly,  the  bulk  of  the  condensed  phase  is  often  assigned  fixed 
properties  and  the  structure  of  the  vapor-liquid(solid)  interface  is  reduced  to  a  surface  E  bounding  the  vapor.  Molecular 
exchanges  through  E  are  modelled  by  the  following  inhomogeneous  linear  boundary  condition  for  the  distribution 
function  /(x,  v|f )  of  molecular  velocities  v  at  a  location  x  on  E: 

(von)/(x,v|f)  =  (von)/e(v)+  [  f?(v,vi)|v,  on|/(x,vi|f)c/3vi,  v°n>0  (1) 

J(\lon)<0 


Eq.(l)  states  that  the  molecular  flux  emerging  from  a  surface  element  of  normal  unit  vector  n  has  two  components. 
The  first  one  consists  of  the  molecules  initially  belonging  to  the  condensed  phase  and  evaporating  into  the  vapor  phase 
with  distribution  function  fe(y),  the  second  one  of  vapor  molecules  whose  initial  velocity  vj  is  instaneously  changed 
to  v  by  the  interaction  with  the  condensed  phase,  described  by  the  scattering  kernel  R(\,\ i).  The  usual  choice  for  the 
distribution  function  /e(v)  of  molecules  evaporating  from  the  condensed  phase  is  the  half  range  Maxwellian 
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being  nw  the  number  density  of  the  saturated  vapor  at  the  temperature  Tw,  whereas  Maxwell’s  gas-surface  scattering 
kernel  is  the  usual  choice  to  describe  molecular  re-emission  from  the  condensed  phase: 
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The  evaporation  coefficient  a,,  (0  <  a,,  <  1)  gives  the  fraction  of  vapor  molecules  impinging  on  the  interface  and 
absorbed.  The  total  fraction  of  impinging  molecules  which  are  instantaneously  re-emitted  is  1  —  ae,  a  being  the 
probability  of  diffuse  re-emission  and  l  —  a  the  probability  of  specular  reflection.  The  assessment  of  the  model 
accuracy  and/or  the  determination  of  the  model  parameters  requires  either  difficult  experiments  or  theoretical  tools 
capable  of  describing  both  phases. 

In  the  present  work  the  evaporation  of  a  fluid  into  vacuum  is  studied  by  the  Enskog- Vlasov  equation[2,  3,  4],  This 
kinetic  equation  describes  a  dense  fluid  whose  molecules  interact  via  an  intermolecular  potential  which  is  splitted  into 
a  hard  core,  treated  as  in  standard  Enskog  equation,  and  an  attractive  tail  that  enters  the  equation  linearly  in  a  mean 
field  term 
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FIGURE  1.  (a)  Density  profiles  and  (b)  Temperature  profiles  for  T  jTc  =  0.596,0.663,0.729,0.795,0.862. 
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being  0(r)  the  Sutherland  potential  and  C£(/i,/i)  the  Enskog  collisional  integral. 

Eq.(4)  provides  a  simplified  description  of  the  microscopic  behavior  of  the  fluid  but  it  has  the  capability  of  handling 
both  the  liquid  and  vapor  phases,  thus  eliminating  the  necessity  of  postulating  ad  hoc  models  for  boundary  conditions 
at  the  vapor-liquid  interface. 


EVAPORATION  INTO  VACUUM:  NUMERICAL  RESULTS  AND  DISCUSSION 

The  behavior  of  the  distribution  function  fe  of  the  molecules  initially  in  the  liquid  phase  and  entering  the  vapor  region 
through  the  vapor-liquid  interface  has  already  been  obtained  in  previous  MD  studies  of  evaporation/condensation 
phenomena[5,  6],  and  deviations  from  the  Maxwellian  behavior  have  been  reported.  However,  in  the  MD  simulations 
cited  above  it  was  necessary  to  adopt  some  criterion  to  distinguish  the  contribution  of  fe  from  the  contribution  of  vapor 
molecules  entering  the  condensed  phase  and  reflected  back  to  the  vapor  region  after  a  short  time  interaction  with  the 
interface  molecules.  In  this  work  the  unwanted  contribution  of  reflected  vapor  molecules  is  minimized  by  placing  a 
perfectly  absorbing  planar  wall,  close  and  parallel  to  the  vapor-liquid  interface^,  8],  The  distance  between  the  wall 
and  the  interface  is  much  smaller  than  the  mean  free  path  in  the  gas  phase  so  that  most  molecules  leave  the  liquid 
surface  and  hit  the  absorbing  wall  without  suffering  intermediate  collision,  but  larger  than  a  few  molecular  diameters, 
in  order  to  avoid  interferencies  with  the  formation  of  the  interface.  A  particle  scheme  [9]  has  been  used  to  calculate 
approximate  solutions  of  Eq(4).  The  spatial  domain  is  a  finite  symmetric  interval  [-L/2.  L/2]  of  the  x  axis  with  two 
perfectly  absorbing  walls  located  at  the  boundary  x  =  ±L/2.  Each  computation  is  started  by  arranging  a  homogeneous 
liquid  slab  at  temperature  7/,  in  the  center  of  the  computational  domain.  The  progressive  cooling  of  the  system  is 
avoided  by  thermostating[10]  at  constant  temperature  Ti  a  thin  strip,  4 a  wide,  in  the  central  part  of  the  slab. 

The  net  evaporation  flow  causes  the  consumption  of  the  liquid  slab  and  the  slow  movement  of  the  vapor-liquid 
interfaces.  Left  to  itself,  the  resulting  flow  would  be  unsteady  and  the  possibility  of  computing  macroscopic  quantities 
by  time  averaging  in  the  poorly  populated  vapor  region  would  be  lost.  However,  it  should  be  observed  that  the  velocity 
of  the  interface  is  much  smaller  than  the  characteristic  vapor  velocity,  at  temperatures  well  below  Tc,  hence  the  flow 
will  be  quasi  steady.  To  limit  the  interface  motion  and  use  time  averaging  of  flow  properties,  the  following  procedure 
has  been  adopted.  After  the  initial  transient  dies  out,  time  averaging  of  system  properties  is  performed  over  a  short 
time  interval.  During  this  time,  the  motion  of  the  interface  is  negligible.  At  the  end  of  this  time  interval,  the  interfaces 
are  reset  in  their  original  position  by  translating  the  fluid  particles  in  each  domain  towards  the  closest  absorbing  wall. 
The  thin  empty  gap,  created  in  the  center  of  the  slab  by  particles  translation,  is  filled  with  the  particles  absorbed  at  the 
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FIGURE  2.  Non-equilibrium  profiles  in  the  vapor-liquid  interface  region;  o,  number  density  na3;  <,  normalized  transversal 
temperature  7j_/7f;  t>,  normalized  longitudinal  temperature  7j|/7),;  A,  normalized  overall  temperature  T / Ti\ o,  normalized  velocity 
uxj mean  field  —1/10 Fxa/(RTi)\  filled  symbols:  DSMC  Bolzmann  equation  solution  with  diffusive  boundary  condition 
at  x  =  xm. 


walls;  their  new  random  velocities  are  generated  by  sampling  a  Maxwellian  with  temperature  7/.  Since  the  gap  width 
equals  the  sum  of  interfaces  displacements  during  the  averaging  time  interval,  the  density  in  the  center  of  the  slab  keeps 
its  equilibrium  constant  value.  At  the  end  of  the  translation  the  fluid  motion  is  again  calculated  and  time  averaging 
continues.  The  results  obtained  by  the  method  outlined  above  have  been  successfully  compared  with  a  simpler,  but 
much  more  time  consuming,  method  in  which  interfaces  are  free  to  move  and  the  unsteady  flow  properties  are  obtained 
from  the  superposition  of  a  few  hundreds  independent  program  runs. 

Simulations  have  been  performed  in  a  relatively  narrow  temperature  range.  The  lower  temperature  limit  is  given  by  the 
necessity  of  keeping  the  liquid  density  within  the  validity  range  of  the  Carnahan-Starling  approximation  and  keeping 
the  vapor  density  high  enough  for  significant  statistics.  The  higher  temperature  limit  is  given  by  the  request  of  having 
a  dilute  vapor  phase.  Density  profiles  are  shown  in  Fig.  1(a).  Each  density  profile  has  a  minimum  in  the  center  of 
the  slab  and  two  local  maxima  in  proximity  of  the  interfaces  because  the  evaporation  cooling  of  the  slab  causes  a 
temperature  decrease  and  a  density  increase  in  the  external  region  of  the  liquid  slab.  The  corresponding  temperature 
profiles  are  shown  in  Fig.  1(b).  The  central  thermostatated  strip  keeps  the  prescribed  temperature  value  7)  .  In  the  liquid 
bulk,  temperature  profiles  exhibit  an  almost  linear  trend  whose  slope  is  more  pronounced  for  the  highest  temperature 
values.  Temperature  profiles  suffer  a  sudden  drop  in  the  interface  regions,  but  the  rate  of  temperature  change  rapidly 
diminishes  in  the  vapor  region.  The  small,  but  evident,  temperature  gradient  in  the  gas  phase  indicates  that  the  effect 
of  collisions  in  the  gas  phase  is  not  completely  negligible. 

The  non-equilibrium  structure  of  the  interface  is  outlined  in  Fig.  2,  where  density,  mean  field  Fx  and  mean  velocity 
ux  are  shown  along  with  the  profiles  of  transversal  (TjJ,  longitudinal  (7y)  and  total  temperature  T  =  (7j_  -F  27’  ) /3. 
The  results  show  that  T±,  ?j|  and  T  follow  the  same  almost  linear  behavior  not  only  in  the  bulk  of  the  liquid  phase 
but  also  in  the  interface  region.  The  three  curves  separate  at  a  position  xs  where  the  density  amounts  to  about  25%  of 
its  maximum  value.  This  separation  position  marks  the  beginning  of  a  transition  layer  which  extends  a  few  molecular 
diameters  into  the  low  density  region.  This  layer  is  characterized  by  the  action  of  the  mean  field  Fx  whose  effects  are 
still  felt  at  a  distance  of  about  3 a  from  the  center  of  the  interface,  defined  as  the  point  xc  where  the  density  takes  the 
value  (n/  +  ng)/ 2.  After  the  separation  position,  T±  is  almost  constant  and  the  structure  of  the  transition  layer  is  better 
shown  by  7’  |  whose  rapid  drop  follows  the  increase  of  ux.  As  shown  by  the  sudden  decrease  in  their  rate  of  change, 
the  velocity  ux  and  the  longitudinal  temperature  7’  complete  their  transition  to  the  dilute  gas  regime  at  a  distance  x,„ 
of  about  4a  from  the  interface  center.  As  discussed  below  in  greater  detail,  this  is  the  position  xm  where  the  matching 
with  Boltzmann  equation  solutions  should  be  done  and  the  boundary  condition  applied. 


499 


FIGURE  3.  (a)  Temperature  at  the  separation  location  as  a  function  of  T Solid  line,  Ts  from  numerical  solutions  of  Eq.(4);  o,  Ts 

from  Eq.(6);  (b)  Evaporation  coefficient  ae  as  a  function  of  Tl-  o,  ae  from  ng[Ti)\  y,  ae  from  ng(Ts)\  A,  oe  from  ng(Ts )  and  back 
scattered  flux  correction. 


The  behavior  of  Ts,  the  temperature  value  at  xs,  as  a  function  of  7/  is  shown  in  Fig. 3(a).  The  separation  temperature 
Ts  increases  more  slowly  than  7/  because  the  evaporation  cooling  becomes  more  and  more  important  as  7/  grows.  In 
the  explored  temperature  range,  the  temperature  difference  7/  —  Ts  exhibits  a  linear  dependence  from  the  mass  flux  Jm . 
The  above  results  suggest  a  simple  but  approximate  method  to  calculate  Ts  by  assuming  that  the  temperature  field  in 
the  liquid  slab  center  up  the  separation  point  can  be  described  by  the  classical  heat  transport  equation.  Accordingly,  at 
the  separation  point  the  following  relationship  holds: 


-X(TS) 


(5) 


where  j(7()  is  the  heat  conductivity  and  A E  is  the  total  energy  per  particle  carried  away  from  the  liquid  phase  by 
evaporation  through  an  imaginary  surface  located  at  xs.  Since  the  flow  is  quasi  stationary  and  the  temperature  field  is 
close  to  a  linear  function  of  x,  the  boundary  condition  (5)  can  be  turned  into  the  following  equation  for  Ts 


Tl-Ts  = 


,  ,  AH  aen,,{Ts) 

X(TS)  gK  J 


(6) 


by  approximating  ^  with  (Ts  —  Tl )  /xs,  A E  with  the  latent  heat  of  vaporization  AH  and  Jm  with  aeng 
heat  conductivity  x(T)  is  not  affected  by  the  attractive  tail[3],  hence  it  can  be  computed  from  its  Enskog  theory 
approximation [1 1],  As  shown  in  Fig. 3(a),  Eq.(6)  slightly  underestimates  7),  but  reproduces  the  correct  behavior  as  Tl 
is  changed. 


EVAPORATION  COEFFICIENTS  AND  DISTRIBUTION  FUNCTIONS 

Simulations  results  show  that,  after  an  initial  transient,  the  particles  flux  at  the  absorbing  wall  reaches  a  constant 
value  Jm  which  can  be  used  to  estimate  the  evaporation  coefficient  ae.  It  should  be  observed  that  the  net  particles  flux 
Jm  at  the  wall  underestimates  the  true  evaporating  flux  Jm  because  a  few  evaporating  particles  are  scattered  back  toward 
the  interface  by  collisions.  However,  in  the  conditions  under  examination,  the  applied  correction  does  not  exceed  2%. 
The  calculation  of  ae  requires  a  reference  Maxwellian  flux  which  can  be  obtained  either  from  Tl  or  from  some  other 
reference  temperature  that  can  be  attributed  to  the  evaporating  liquid  surface.  It  is  not  clear  how  a  specific  temperature 
value  can  be  assigned  to  the  liquid  surface  in  the  present  non-equilibrium  conditions,  however  the  temperature  Ts  is 
also  a  reasonable  choice  as  a  reference  temperature,  since  it  takes  into  account  the  effects  of  evaporation  cooling  and  it 
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FIGURE  4.  (a)  Reduced  distribution  functions  of  particles  absorbed  at  walls;  □,  fx{vx);  A,  fy ( vy ) ;  \7,  fz(vz)\  dashed  line, 

half  and  full  range  Maxwellian  at  temperature  0.96 77;  (b)  Reduced  distribution  function  fx(vx)  at  x  =  x,n  as  a  function  of  77; 
o,  77/77  =  0.596;  □,  77/77  =  0.0.663;  0,77/7/  =  0.862;  dashed  line,  half  range  Maxwellian  with  unit  temperature.  Velocities 
normalized  to  ^ /RTref ,  Tref  =  Tj_(xm). 


is  possible  to  estimate  it  from  a  simple  macroscopic  model  of  the  liquid  phase.  Accordingly,  values  of  the  evaporation 
coefficient  oe  have  been  computed  from  the  corrected  flux  Jm  using  two  distinct  reference  temperatures 


®etlg(T) 


(  T  =  Tl 
\t  =  Ts 


(7) 


As  shown  in  Fig3(b),  the  computation  of  ae  from  7/  leads  to  a  strong  temperature  dependence  of  the  evaporation 
coefficient  because  the  evaporation  cooling  of  the  slab  strongly  reduces  the  evaporation  flux.  However,  if  Ts  is  used, 
then  the  temperature  variation  of  the  evaporation  coefficient  is  very  much  reduced  and  values  of  ae  around  0.9  are 
obtained  in  all  of  the  cases  examined,  as  also  observed  in  Ref.  [8].  As  briefly  mentioned  above,  if  collisions  in  the 
vapor  phase  could  be  neglected,  then  fe  could  be  reconstructed  from  the  particles  absorbed  at  the  walls.  Fig. 4(a)  shows 
the  reduced  distribution  functions  of  the  velocity  component  vx,  normal  to  the  walls  and  of  the  velocity  components  vy 
and  vz,  parallel  to  the  walls,  denoted  by  fx(vx),  fy(vy )  and  fz(vz),  respectively.  The  numerical  distribution  functions  are 
normalized  to  unity  and  compared  with  a  Maxwellian  distribution  function  with  zero  bulk  velocity  and  a  temperature 
obtained  from  the  vy  and  vz  velocity  components.  It  is  clear  that  the  distribution  functions  fy(vy)  and  fz(vz)  are 
very  well  approximated  by  a  Maxwellian,  whereas  the  distribution  function  fx(vx)  shows  a  deviation  from  the  half¬ 
range  Maxwellian  that  one  would  observe  if  molecules  were  emitted  from  the  liquid  surface  according  to  Eq.(2).  The 
deviation  occurs  in  the  low  velocity  region  and  it  is  due  to  collisions  in  the  vapor  phase  which  preferentially  remove 
slow  particles  that  take  a  long  time  to  reach  the  absorbing  wall.  The  effect  of  temperature  on  the  shape  of  fx{vx)  at 
xm  is  described  by  Fig. 4(b).  The  results  show  that  for  small  velocity  values,  a  deviation  from  the  local  half-range 
Maxwellian  is  observed;  the  difference  grows  as  7/  increases.  It  should  be  observed  that  the  velocity  distribution 
functions  described  above  do  not  take  into  account  the  history  of  individual  molecules.  As  a  matter  of  fact,  the  positive 
velocity  range  of  fx(vx )  at  a  certain  location  in  the  interface  region  will  mostly  contain  the  contribution  of  molecules 
originated  in  the  liquid  region  and  travelling  toward  the  vapor  region.  However,  it  will  also  contain  a  contribution  from 
vapor  molecules  which  have  been  reflected  back  to  the  vapor  phase.  The  effects  of  this  second  molecular  group  is 
small,  but  the  deviation  from  the  Maxwellian  behavior  is  also  small  and  it  is  necessary  to  check  whether  the  effect 
is  due  to  the  reflected  component.  To  distinguish  the  two  molecular  groups  a  “color”  variable  has  been  added  to 
each  simulation  particle.  Particles  in  the  strip  |x|  <  xm  initially  get  the  same  “color”.  During  a  simulation,  the  color 
variable  of  a  particle  is  changed  only  once  on  entering  the  vapor  region  from  the  interface-vapor  boundaries  at  ±xm. 
As  expected,  the  negative  tail  of  the  distribution  function  of  molecules  which  have  never  crossed  the  interface-vapor 
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boundaries  is  less  populated  but  the  shape  in  the  positive  velocity  range  is  the  same  as  the  full  distribution  fx(vx), 
indicating  that  the  deviation  from  Maxwellian  behavior  is  exhibited  by  the  evaporating  component  and  it  is  not  due  to 
the  small  contribution  of  the  reflected  component. 

The  results  described  above  indicate  that,  at  the  boundary  between  the  vapor  phase  and  the  interface,  the  distribution 
function  of  evaporating  molecules  fx(vx)  deviates  to  some  extent  from  Eq.(2).  Since  the  deviation  is  rather  limited,  it  is 
quite  natural  to  ask  whether  the  vapor  motion  between  the  interfaces  at  ±xm  and  the  absorbing  walls  can  be  described 
by  the  Boltzmann  equation  with  boundary  conditions  of  the  kind  given  by  Eqs.(2)  and  (3).  To  answer  the  question  a 
few  DSMC  simulations  have  been  performed  to  study  the  motion  of  a  dilute  hard  sphere  vapor  in  the  gap  between 
xm  and  the  absorbing  wall.  The  gap  width  is  L/  —xm.  The  fluid  at  the  left  hand  side  of  the  interface  position  has  been 
replaced  by  a  purely  diffusive  boundary  condition  (a  =  1)  with  ae  =  0.9.  The  saturated  vapor  density  nw  has  been 
taken  equal  to  ng(Ts )  and  the  wall  temperature  Tw  has  been  taken  initially  equal  to  T±(x,„)  and  then  adjusted  to  fit  the 
temperature  profiles  in  the  gas  phase  obtained  by  the  Enskog-Vlasov  equation  solutions.  As  is  clear,  the  procedure  to 
match  Boltzmann  equation  to  Enskog-Vlasov  equation  solutions  is  not  rigorous  at  all.  However,  the  comparison  can 
give  an  indication  of  the  accuracy  of  the  boundary  condition  models.  Fig. 2  shows  the  comparison  of  temperatures  and 
velocity  profiles  in  the  vapor  phase  for  Ti/Tc  =  0.596  and  L/  8a.  This  case  corresponds  to  the  lowest  temperature 
value  where  good  agreement  has  been  found  between  fx{vx)  and  a  half-range  Maxwellian  for  vx  >  0.  The  temperature 
profiles  obtained  by  the  numerical  solution  of  the  Boltzmann  equation  are  in  excellent  agreement  with  the  Enskog- 
Vlasov  profiles.  The  velocity  profile  computed  from  Boltzmann  equation  is  slightly  above  the  Enskog-Vlasov  profile. 
A  closer  inspection  of  the  solutions  data  shows  that  the  discrepancy  is  due  to  the  negative  tail  of  fx  which,  in  the 
Boltzmann  equation  solution,  is  less  populated  than  its  Enskog-Vlasov  counterpart.  The  comparison  of  density  profiles 
shows  that  nw  =  ng(Ts)  is  a  reasonable  choice  in  this  particular  case.  For  higher  temperature,  for  which  the  largest 
deviation  from  the  half-range  Maxwellian  has  been  observed,  the  agreement  is  obviously  less  good. 


CONCLUSIONS 

The  present  work  is  an  attempt  to  give  a  unified  description  of  a  two-phase  system  through  an  approximate  kinetic 
equation.  The  evaporation  of  a  liquid  slab  into  near  vacuum  is  investigated  to  provide  information  about  the  distribution 
function  of  molecules  emitted  from  the  liquid  phase.  The  non-equilibrium  structure  of  the  vapor-liquid  interface  for 
various  temperature  values  has  been  obtained.  In  particular,  the  profiles  of  macroscopic  quantities  in  the  transition 
layer  bridging  the  liquid  to  the  vapor  phase  has  been  studied.  The  examination  of  reduced  distribution  functions  shows 
that  the  velocity  component  parallel  to  evaporating  surface  is  very  well  approximated  by  a  full  range  Maxwellian.  The 
distribution  function  of  the  normal  velocity  component  of  outgoing  molecules  is  close  to  a  half-range  Maxwellian  at  the 
lower  values  of  the  investigated  temperature  range,  but  deviates  from  a  half-range  Maxwellian  at  higher  temperature. 
These  results  agree  with  the  recent  MD  studies  of  the  same  problem[12],  where  a  similar  behavior  of  the  distribution 
function  has  been  found  in  the  evaporation  of  the  Lennard-Jones  fluid.  Therefore,  the  Enskog-Vlasov  model  seems  to 
provide  a  sufficiently  realistic  description  of  two-phase  flows,  in  spite  of  its  approximate  nature. 
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